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Abstract. We show that the full dynamical freedom of the well known Szekeres models al¬ 
lows for the description of elaborated 3-dimensional networks of cold dark matter structures 
(over-densities and/or density voids) undergoing “pancake” collapse. By reducing Einstein’s 
field equations to a set of evolution equations, which themselves reduce in the linear limit 
to evolution equations for linear perturbations, we determine the dynamics of such struc¬ 
tures, with the spatial comoving location of each structure uniquely specified by standard 
early Universe initial conditions. By means of a representative example we examine in de¬ 
tail the density contrast, the Hubble flow and peculiar velocities of structures that evolved, 
from linear initial data at the last scattering surface, to fully non-linear 10-20 Mpc scale 
configurations today. To motivate further research, we provide a qualitative discussion on 
the connection of Szekeres models with linear perturbations and the pancake collapse of the 
Zeldovich approximation. This type of structure modelling provides a coarse grained - but 
fully relativistic non-linear and non-perturbative - description of evolving large scale cosmic 
structures before their virialisation, and as such it has an enormous potential for applications 
in cosmological research. 
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1 Introduction 

Considering numerical solutions of Einstein’s equations applied to Cosmology is evidently an 
urgent task, as it is widely assumed to be practically impossible to model minimally realistic 
cosmic structures (even at a coarse grained level) by means of exact solutions of Einstein’s 
equations. Since cosmological application of numerical General Relativity is still in its early 
stage of development [1—3], most cosmological applications that require a non-perturbative 
relativistic approach still rely on the highly idealised class of spherically symmetric Lemaitre- 
Tolman-Bondi (LTB) dust models [4, 5], which can only describe the evolution of a single 
spherical dust structure embedded in an FLRW background. 

The non-spherical Szekeres dust models (see details of their classification in [4, 5]) are 
a well known generalisation of LTB models. Although it is wholly unreasonable to expect 
of these models (themselves an exact solution Einstein’s equations) to provide the level of 
“realism” expected from the (yet to develop) numerical solutions, Szekeres models are still 
useful and have been applied to address various cosmological and observational issues [6-22], 
However, in practically all this literature the authors consider models that only describe the 
evolution of two structures (an over-density next to a density void) in simple axial dipolar 
arrays. In a recent article [23] we showed that the full dynamical freedom of the models allows 
for the description of far more general configurations, namely: elaborated networks of multiple 
evolving cosmic structures (over-densities and density voids defined by 3-dimensional maxima 
and minima of the density), whose spatial (radial and angular) location at all t can be a priori 
specified by suitable initial conditions. However, [23] was essentially a theoretical study 
concerned with the existence conditions of the extrema (maxima, minima and saddle points) 
of all Szekeres scalars (not just the density) of generic models (ever expading, collapsing, 
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with zero and nonzero A). While various relevant technical issues, such as avoidance of shell 
crossings, were extensively discussed, this study only provided a simple qualitative example 
(see its figures 8 and 9) of the density contrast of a multi-structure Szekeres configuration 
at an initial time, without studying its actual evolution and without illustrating the shape of 
other scalars (for example the Hubble scalar). 

In the present article we aim at extending and enhancing the work in [23], specifically 
by undertaking the following tasks: (i) Modelling more general structure networks. 
The structures in [23] only admitted an over-density or a density void in each radial shell 
zone. More general networks of structures are now obtained, admitting an arbitrary number 
of structures in assorted angular locations in each radial zone. With this improvement we can 
now attempt to implement a coarse grained modelling of large scale cosmography, either from 
observed and/or reconstructed studies [24, 25] or from numerical simulations [26]. (ii) Real¬ 
istic evolution. Assuming a ACDM background consistent with observations, we examine 
the numerical evolution of the above mentioned improved networks, from linear perturbations 
at the last scattering surface into an ~ 80 Mpc region containing 10-20 Mpc sized structures 
in a non-linear regime at present cosmic time, (iii) Expansion, collapse and peculiar 
velocities. We obtain and depict the present day anisotropic and inhomogeneous Hubble 
flow that shows the structures undergoing “pancake” collapse at rates appropriate to their 
length scales. The peculiar velocities of the structures are examined, providing a qualitative 
comparison with velocities reported in the existing literature, (iv) Theoretical issues. We 
show that the pancake collapsing Szekeres over-densities provide an exact relativistic ana¬ 
logue of the Newtonian Zeldovich Approximation. We also comment on the correspondence 
with linear perturbations of dust sources (in the isochronous gauge) and on the back-reaction 
issue. 

The section by section contents are as follows. In section 2 we derive a set of evolution 
equations (equivalent to Einstein’s held equation) that fully determine the dynamics of Szek¬ 
eres models. A procedure to construct Szekeres configurations describing elaborated networks 
of structures (over-densities and density voids) is presented in section 3. In section 4 we build 
up a representative numerical example consisting of a central spheroidal void surrounded by 
multiple over-densities undergoing pancake collapse. In section 5 we examine in detail (i) 
the density contrast, the Hubble scalar and eigenvalues of the expansion tensor and radial 
peculiar velocities associated with the example of section 4. The connection with the Zel¬ 
dovich approximation and dust linear perturbations are discussed in section 6. Conclusions 
and guidelines for further research and applications are stated in section 7. 


2 The dynamics of Szekeres models 


Quasi-spherical Szekeres models of class I in “stereographic” spherical coordinates (r, 9, <fi) 
[4, 6, 7] are described by the following non-diagonal metric 1 



(T-W) 2 sin 4 9 

1 - [tCqjinir 2 + (1 + cos/9) 2 


w 1 - 2 1 + i os - Z W 


sin 


9re — 


1 + cos 9 
.22 


(W - Z ), g rcf> = - 


a 2 r sin 2 6 
1 + cos 9 


W, 


gee = a~r~, = a 2 r 2 sin 2 9, 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 


1 A simpler diagonal metric (with several variations) is used in most of the existing literature. The trans¬ 

formations relating (2.1)—(2.3) to this metric are given in Appendix A of [23]. 
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where a = a(t,r), T = 1 + ra'/a, with a’ = da/dr, [/C 9 ]; n i = [JC 9 ]mi(0 (see (2.13)), and the 
Szekeres dipole W is given by 


W = —X sin 6 cos (f> — Y sin 6 sin 4> — Z cos 9, (2.4) 

where X = X(r), Y = Y(r), Z = Z(r) are the dipole free functions and W = VX 2 + Y 2 + Z 2 
is the dipole magnitude. It is straightforward to see from (2.1)—(2.3) that the surfaces of 
constant t and r are 2-spheres (with surface area 47ra 2 ?’ 2 ) that are non-concentric about the 
origin worldline [4, 23]. By setting X = Y = Z = 0 =+ W = 0 we obtain a generic LTB 
“seed model” as the unique spherically symmetric sub-case. 

The models are fully characterised by their covariant fluid flow scalars: the density p, 
the Hubble scalar % = 0/3 and the spatial curvature JC = (1/6)^TZ (with © = V a u a and 
ZPfZ the Ricci scalar of the hypersurfaces of constant t). Considering a nonzero cosmological 
constant to accommodate a ACDM background, the dynamics of the models becomes fully 
determined by the numerical solutions of the following evolution equations derived in [17]: 2 


Pq 

— 3p q 7~Lq, 

(2.5) 

ii q 

- , 2 4tt 8tt 

(2.6) 

A (p) 

= -3(1 +A (p) )D (w) 

(2.7) 

d (h) 

= (-27A, + 3D^>) 

(2.8) 

a 

— CL 'hiq , 

(2.9) 

Q 

= 3£D (W) , £ = r_W , 

y y i _ w 

(2.10) 


subject to the algebraic constraints: 



T~L 2 q — — [p q + A] — ICq, 

(2.11) 


2U q T> {n) = — - D^), 

3 

(2.12) 

where the “q-scalars” 
covariant scalars) are 

A q and their exact fluctuations [17] (which determine the standard 

given by 


„ [P?]ini ^ [^Cg]ini qj 0, 

Pq — q 5 — 9 5 T~L-q — » 

a z a 

(2.13) 


°'"- A A --3(r -wy A- p ,u,k, 

A (p) A(P) P~Pq r Pq/P 9 

Pq Pq 3(T- W)’ 

(2.14) 

(2.15) 


with the subindex j n ; denoting henceforth evaluation at an arbitrary time slice t = t m \. 

Besides the cosmological constant A, the initial conditions to integrate the system (2.5)- 
(2.8) are the following five free parameters: the two “radial” initial functions common to 

2 In practically all the existing literature the dynamics of Szekeres models is studied in terms of the integral 
solutions (analytic or numerical) of the Friedman-like equation (2.11). See for example the comprehensive 
work in [18], which only considered models with A = 0. We believe that the evolution equations (2.5)-(2.10) 
provide a much more efficient framework for numerical work, specially for models with A > 0. 
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LTB models, [/? g ]mi, [/C„]i n i, and the “angular dipole” initial functions X, Y, Z. Initial values 
of Hg, AM, D^, follow from (2.11)—(2.12), the radial coordinate is chosen so that 

dini = Tin; = 1, while the Big Bang time and its gradient follow from the choice of [/9q]i n i, [fcq]mi 
(see [17, 23]). 

The q-scalars and their fluctuations are coordinate independent objects that are directly 
related to curvature and kinematic scalars [17, 23]. As we show in section 6.2, they reduce in 
the linear limit to standard variables of cosmological dust perturbations in the synchronous 
gauge (see the LTB case in [27]). 

3 Networks of over—densities and density voids 

Over-densities and density voids can be defined as regions surrounding the spatial maxima 
and minima of the matter density. The coordinate location of the density extrema (as that 
of all other scalars A) follows from the condition A' = A g = A ^ = 0, whose solutions are, at 
each constant t, 

r = r e ±, 9 = 6±(r e ±), 0 = (f>± (r e ±) (3.1) 

where 9±{r) and <j>±(r) follow from the solutions of the subset Ag = A ^ = 0 (the “angular 
extrema” [23]), 


= arctan 



0+ = 7T + 


= arccos 


\JX 2 + Y 2 + Z 2 


9 + = 7T — 6 -, 


(3.2) 

(3.3) 


which for every fixed r defines a precise angular direction and also, for varying r, the two 
“curves of angular extrema” B±(r) = [r , d±(r), (j>±(r)] in all time slices. To find the radial 
location r = r e ± of the extrema we need to solve the “radial” conditions A'±{t,r) = 0, where 
the subindex ± denotes evaluation along the curves B±[r) (notice that W± = ±W). 

For each solution of A!_ j_ = 0 at arbitrary t there will be an extremum of the scalar A at 
angular coordinates (3.2)-(3.3). As shown in [23], a sufficient condition for the existence at 
all t (pending shell crossings) of an arbitrary number of such solutions (and thus an arbitrary 
number of spatial extrema of all Szekeres scalars) is furnished by assuming compatibility with 
Periodic Local Homogeneity (PLH), defined by the vanishing for all t of the shear (a£) and 
electric Weyl (E£) tensors along a sequence of comoving 2-spheres (comoving homogeneity 
spheres) [23]. Since (2.12) and (2.14) are preserved by the time evolution, models compatible 
with PLH are specified by the following initial conditions: 3 


=► D 1 ( n f(r:,0,0) = O ^ [A' q } ini (rt) = 0, 

=> Ani (rl,0,4>) = [A q \i ni (K), (3.4) 


which imply that [a£]* = [-£)(]* = 0 holds for all t, with the subindex * denoting evaluation 
at the sequence of n nonzero values of the radial coordinate r\. i = 1 ,..,n that mark the 
comoving homogeneity spheres. 

3 It is important to emphasise that PLH is a sufficient (but not necessary) condition for the existence of 
spatial maxima and minima of Szekeres scalars. These extrema can also arise without assuming PLH by 
means of “simulated shell crossings” induced by the dipole parameters X, Y, Z for arbitrary choices of initial 
value functions [pq]i n i, [^CJini, [Hq] ini- See detail in [23]. 
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Depending on the number of values r\ in the comoving shells defined by the intervals 
A* = r* -1 < r < r\. initial conditions (3.4) lead to 2n + 1 spatial extrema distributed in the 
A* as follows (see comprehensive discussion in [23]): 

• If regularity conditions hold [4, 23] the origin r = 0 is always a spatial maximum or a 
spatial minimum (depending on the sign of A" at r = 0). 

• There are n maxima or minima along (this depends on the profile of A± and the sign 
of D (i4 >, see figure 1), while the n extrema along B-(r) are necessarily spatial saddles. 

However, as we show in this article, the models admit more general configurations with an 
arbitrary number of spatial maxima or minima at each radial interval A*. Further restrictions 
on the functions [H g ]; n i satisfying conditions (3.4) may be necessary to ensure the absence of 
shell crossings and the preservation of the initial concavity for all the evolution (see sections 
IX and X of [23]). 

It important to emphasise that the comoving shells marked by the radial interval A* 
are not “FLRW regions”, since = E? = 0 only holds in the boundaries r = r\ of these 
shells (the “comoving homogeneity spheres”). Since Szekeres models can always be matched 
(along comoving 2-spheres) to regions of dust FLRW models [12], any one of the comoving 
shells marked by some A* can be replaced by a FLRW shell region matched to contiguous 
shells (so that = E£ = 0 holds in the whole interval A*), but this type of matching is too 
restrictive and thus has not been considered. However, even if we had introduced such FLRW 
shell regions, the resulting configuration would bear no resemblance to Swiss Cheese models, 
as in the latter the FLRW “cheese” is not distributed in spherical shells but is surrounding a 
collection of inhomogeneous dust vacuoles (the “holes”). Further discussion on these issues is 
given in the conclusions section. 

4 Numerical example of multiple structures. 

We illustrate the set-up of a Szekeres model describing multiple evolving structures through a 
simple idealised numerical example (more elaborated examples can easily be obtained along 
these lines). For this purpose, we assume an asymptotic ACDM background characterised by 
the present day parameters from Planck 2013 [28]: = 0.32 (includes baryons),HQ = 0.68 

and Hq = 68km/s Mpc (over-bar denotes background ACDM variables). We consider the 
model evolution (governed by (2.5)-(2.12)) from linear initial conditions (see Table 1) at the 
last scattering surface (LS) t m \ = t LS ~ 4 x 10 5 ys that also comply with PLH conditions (3.4). 


4.1 Radial location of the spatial density maxima 

We consider the dimensionless initial density and spatial curvature q-scalars (2.13) as 

[/A?]ls(x) — 3 g 2 ’ D%]ls(x) — , (4.1) 

where H hS ~ 2/(3t LS ) and y = r/r s with r s = 0.0025 Mpc, which fixes the unit comoving 
length scale well within the comoving horizon scale at LS. An initial density minimum at 
r = 0 follows from the condition [P(j]l S ( 0) > 0. For a sequence of density maxima inside 
intervals A* along the curve B + (r) (as in figure 8 of [23]), the profiles of [/r g ] L s and [n q ] LS < 0 
for all y must correspond to non-decreasing functions complying with (3.4) and with a ACDM 
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Ai = 0 < x < X* 

Mls — 2o,i 

x), Mls = Ml (x) 

A0i k 

A 

y 

0n < 0 < 012 

— COS 01101101 

-sin 0110 H 01 

012 < 0 < 2vr + '0H 

COS 01201201 

-sin 0 i 20 i 20 i 

A* = xi < x < xi 

Mls = Qi, 2 (x), Mls = Pi, 2 (x) 

A02fc 

X 

Y 

021 < 0 < 022 

- COS 02102102 

- Sin 02102102 

022 < 0 < 2n + 021 

- COS 02202202 

-sin 02202202 

A 3 = Y 2 < Y < Y 3 

— * A* ^ A ^ A* 

Mls = Q2,3 

X), Mls = ^2,3 (X) 

A03A; 

X 

Y 

031 < 0 < 0 32 

- COS 03103103 

-sin 03103103 

032 < 0 < 033 

- COS 03203203 

- sin 03203203 

033 < 0 < 27T + 03i 

- COS 03303303 

- sin 03303303 

A* = X* < X < X* 

M]ls — 03,4 

X)> Mls = 703,4 (X) 

A04fc 

X 

Y 

041 < 0 < 042 

- COS 04104104 

— sin 04104104 

042 < 0 < 043 

- COS 04204203 

— sin 04204204 

043 < 0 < 27T + 04i 

- COS 04304304 

- sin 04304304 


Table 1. Initial conditions at last scattering. The table displays the piecewise defini¬ 
tion of the free functions [^ 9 ] LS , [/%]ls ; X and Y needed to integrate the system (2.5)-(2.12) 
(as we have assumed Z = 0 and a = T = 1 at i = t* = f LS ). The fifth order poly¬ 
nomials Qj-i,k and Vj-i t k are defined by conditions (4.2). The normalised coordinates of 
the comoving homogeneity spheres are [xi,X*,xi,X*] = [3.47826,6.06906,7.99883,9.43622], 
with x* = 0- The azimuthal angular location of the maxima are given by 

[0n, 012,021,022, 031,032,033, 041,042,043] = [0, 57r/4, 37t/4, 77t/4, 0, 37t/3, 47t/3, 7t/3, 7T, 57t/3]- The 
boundaries of the azimuthal partitions are [V’n, 0i2, 021, 022, ^31 , 032, 033, 04i, 04i, 042, 043] = 

[1.96466,5.10625, 7t/4, 57r/4, 27t/6, 7r, 57t/3, 0,27r/3, 47t/3]. The amplitude constants are 

C 11 = 0.885, Ci 2 = 0.89, C 21 = C 22 = C 23 = 0.823, C 31 = C 32 = C 33 = 0.844 and 
C4i = C 42 = C 43 = 0.8592. 


background, whose profiles are depicted by figure 1 (see also panel (a) of figure 4 in [23]). 
We specify [Mls and [Mls as piecewise functions defined at each interval A® for a sequence 
of four intervals x* with x* = 0 (see Table 1), where Q*-i,i(x) and Pi— i,*(x) are fifth order 
polynomials whose six coefficients are determined (at each A®) by two boundary conditions 
and four conditions to fulfill smoothness (of the metric, the covariant scalars and their first 
derivatives) at the dimensionless comoving radii xt : 


Qi-iAxt 1 ) 
Pi-i,i (xi” 1 ) 


m(xl 1 ), 

Mxr 1 ), 


Qi-i,i(xl) = Mxl), 

Pi-i,j (xi) = Hx i), 


O' 1 ■ = O" , . = 0 
V 1 , . = V" , . = 0 

' i—\,i 7 i— l,z ’ 


at X = xl \xl, 
at X = xlMxl, 


(4.2) 
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Figure 1. Profiles of the “radial” initial value functions [/r 9 ] LS and [«-, y ] LS . These functions, 
defined by (4.2) in terms of fifth order polynomials at each interval A*, are depicted as solid black 
curves (only two intervals A* are displayed). The dashed curves are the auxiliary functions m(x) and 
k(x) defined in the text and the blue curves are the initial functions // LS and « LS evaluated along the 
curves £>+(r), which are a collection of piecewise continuous line segments (see figures 2, 3 and 4). 
The maxima of these curves provide the radial coordinate of all maxima at every A*. 




x [Mpc] 


x [Mpc] 


Figure 2. Density contrast. Equatorial projection of the density contrast 6 at initial time f LS (left 
panel) and at present cosmic time t 0 (right panel). The vertical and horizontal scales correspond to 
[x,y\ = [ Rcoscj ), f?sin^>], with R = ar (notice that a = 1 at t = t bs ). Solid radial line segments are 
the curves B+(r). Dashed line segments mark the boundaries of the angular partitions A (j>ik- Dashed 
circles denote the radial comoving values xl that mark the comoving homogeneity spheres. Notice 
that in most of the volume the density contrast is roughly the background value S « 0 (p ~ p) . 

with the auxiliary functions m(x) = 0.5 —0.3/(1+ x 3 ) and k(x) = —0.0014/(l + x 7 ^ 5 ) (dashed 
black curves in figure 1). The normalised radial coordinates Xe+ °f the initial density maxima 
in each A* are the maxima of the curve [/r + ] LS in the left panel of figure 1. To ensure a 
ACDM background at f LS we set [/r g ] LS = m(x) and [«: g ] LS = k(x) f° r X > Xi> leading to 
2 [/Aj] ls —> = 1 and [K g ] LS —> = 0 as r -too (in our setup 3 ~ 10^ 9 ). The functions 

[/r 9 ] LS and [Kq] L s comply with the following desirable properties: (i) they are consistent with 
the conditions to avoid shell crossings for the time range t LS < t < to] (ii) they produce a 
non-simultaneous Big Bang time (t' bb ^ 0), but with negligible differences (of order ~ 10 3 ys) 
in the cosmic age for all observers at to, and (iii) the concavity of the central void and the 
density maxima is preserved for all t > f LS (these technical issues are discussed in detail in 
[23]). 
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4.2 Angular location of the spatial density maxima 

Configurations containing n density maxima (one in each interval A® as in figure 8 of [23]) fol¬ 
low if we assign to each maxima the angular coordinates (0j, <ff) by prescribing at each A* the 
dipole parameters in a piecewise manner: X = — cos </>jSin OiCi fi, Y = — sin 0jSin 9 % Q f \, Z = 
— cos Qi (i fi with the constants 0 < Q < 1 controlling the density contrast amplitude. The n 
functions fi(x) are thus given by 


fi(x) 


sm 


(x-xl )rr 


L Xi - X 


i— 1 


(4.3) 


and satisfy the following boundary and smoothness conditions: /i(x* _1 ) = fi{x*) = 0, /j'(x* _1 ) = 
f[ (x*) = 0 and r *) = = 0. Configurations that are more general than those exam¬ 

ined in [23], admitting several maxima in assorted angles (9 lk , faff) within every A®, follow by 
defining X, Y, Z (at each A®) as piecewise functions (see Table 1): X = — cos 4>i k sin 9i k Cik fi, Y = 
- sin (f) ik sm 9 ik ( ik fa, Z = - cos 9 ik ( ik fi on a partition A<j> ik = ip ik < 4> ik < if ik+1 (with 
k = 1, ,.,p) of angular domains (at each A*) separated by fixed azimuthal angles ipi k whose 
value is chosen to fulfil smoothness conditions between each angular domain b The radial 
coordinate location is the same for all maxima in the same interval A® and the constants 
0 < Q k < 1 define the density contrast amplitude of the maxima. 

For illustrative purposes we select Z = 0, so that (from (3.2)-(3.3)) we have 9i k = 7 t /2 
and thus all spatial density maxima (whose existence is guaranteed by the choice of [p<j]ls 
and [k 9 ] ls ) are located in the equatorial plane 9± = n/2 (the more general case Z / 0 is 
analogous). We select the parameters X and Y in the piecewise manner explained above and 
shown explicitly in Table 1. Since f^ixi) = /l(x*) = 0, we can choose X = Y = 0 (and thus 
W = 0) for y > x*, from this radius the configuration becomes spherically symmetric and 
convergent at all f to a ACDM background [23]. 


5 Discussion. 

By integrating the system (2.5)-(2.12) for the initial conditions specified in Table 1, we can 
examine relevant dynamical quantities that characterise these multi-structure configurations. 

It is important to mention that these configurations are not spherically symmetric, hence 
the apparent rotational symmetry in the resulting graphics displayed in figures 2, 3 and 4 is 
merely an effect arising from employing spherical coordinates. This effect disappears when 
the figures are plotted in terms of proper radial distance (see figure 1 of [23]) or luminosity 
distance. 

5.1 Density contrast. 

We obtain the density contrast 5 = (fi—jYj/p,, where /x = {Att p)/(2>1~L 2 ) and jl(t) = {Airp)/(2>1~L 2 ) = 
Cl m (t)/2, with the Szekeres density obtained from p = p q [ 1 +A( p )] and p(t), il(t) are the den¬ 
sity and Hubble scalar of the ACDM background. Figure 2 displays the level curves of the 

4 The angular partition at each shell A* is equivalent to matching several regions of separate Szekeres 
models with different dipole parameters along common surfaces marked by constant <j>, hence parametrised 
by ( t,r,8). Under certain algebraic restrictions on C ik,4 > ik,Qik the metric and its derivatives tangent to the 
matching surfaces are continuous, hence this matching can be smooth (Darmois conditions hold) even if the 
derivatives with respect to (j> are discontinuous. For the case Z = 0 considered in the numerical example 
Darmois conditions hold without further restrictions. 
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x [Mpc] 



0.009 0.01 0.011 0.012 0.013 0.014 0.015 

r 


Figure 3. The expansion scalar and the eigenvalues of the expansion tensor. The figure 
depicts the equatorial projection of the Hubble scalar H (left panel) evaluated at to (in units km/(s 
Mpc)) and the eigenvalues H^, H 42 ) = H 4 ' 3 ^ of the expansion tensor (right panel) given by (5.1) 
evaluated along a radial ray (curve with t, 9, <f> constant), for t = to/2, to and with 9 = 9i\ = 7t/2, </> = 
</> 2 i = 37 t/ 4 marking the angular coordinates of one of the over-densities listed in see Table 1 (similar 
curves result for all over-densities). Notice that TL as well as H 42 * = H 43 ' are everywhere positive, 
taking (as expected) larger values in the void region than in the over-densities. However, H 414 becomes 
negative at t = to along the over-densities, thus indicating that the latter are undergoing a “pancake” 
collapse. 



Figure 4. Peculiar velocities. Equatorial projection of the radial peculiar velocity field (in km/s) 
at the present cosmic time to, with respect to the background identified as the CMB frame (left panel) 
and with respect to the observer at the centre of the void (right panel). 


equatorial projection of 6, as functions of the area distance R = a(t, r)r, at the initial time 
t = t LS (left panel) and at present time to = 13.7 Gys (right panel). Both panels reveal 
in most of the spherical volume a slight under-density 6 < 0 with near background density 
(6 ~ 0), together with well defined and localised structures: a spheroidal density void around 
the origin (blue shading) surrounded by ten elongated over-densities (red/yellow shading), 
each one around a local density maximum located in the azimuthal angles given in Table 1, in 
the curves £>+(r) (solid line segments) in each one of the angular partitions in each of the four 
intervals A*. As shown in [23], the over-densities have a pancake shape 3-dimensional mor¬ 
phology. Notice that the initial multi-structure shape is preserved in time, but it expanded 
from ~0.08 Mpc at t LS to about ~ 80 Mpc at to, the negative amplitude of the density contrast 
of the void increased three orders of magnitude from <5 LS ~ —2 X 10” 4 at t LS to 5q ~ —0.4 
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at t = to, while the density contrast amplitudes of the over-densities evolved from the linear 
value h LS ~ 3.7 x 10 -4 to fully non-linear values <5o ~ 1.5 — 3, as expected for ~ 10 — 20 Mpc 
sized structures. 

5.2 Expansion and collapse of Szekeres structures. 

The criterion for local collapse in inhomogeneous models follows from the signs of the eigen¬ 
values ((A) = 1, 2,3) of the expansion tensor H^ = Hh < / J -\- a£ (see equations (5)-(6) of 

[29]). Since the shear tensor is traceless, the Hubble scalar H = (1/3)(H (1 ) ± H (2 ^ + H (3 )) 
is the simple arithmetic average expansion. The isotropic spheroidal collapse (three negative 
eigenvalues) implies H < 0, but “pancake” collapse (one negative eigenvalue) and “filamentary” 
collapse (two negative eigenvalues) can occur with H > 0 (overall average expansion). 

For Szekeres models the eigenvalues of H 1 / take the form 5 

H (1) =77 + 2D (w) ='H (? + 3D ( ' h) , H (2) = H (3) = H - D (w) = H q , (5.1) 

where we remind the reader that H q = a/a and = (1/3 )Q/Q follow directly from 

(2.9)-(2.10). It is straightforward to obtain the eigenvalues H ( ^ given by (5.1), evaluated at 
t = to, from the numerical solution of the system (2.5)-(2.12). We the depict in Figure 3 the 
equatorial projection of the Hubble scalar H plotted in terms of the area distance Ro = ao r 
(left panel), as well as the eigenvalues H q 1 "* and Hq 2) = Hq 3 '* evaluated for fixed t = to/2 and 
t = to along a radial ray intersecting one of the over-densities of figure 2. All these quantities 
are given in units of km/(s Mpc). The left panel of figure 3 reveals that most of the displayed 
volume expands at a slightly larger but almost background value Ho ~ Ho = 68km/(s Mpc). 
Both panels also reveal a strong anti-correlation between the Hubble flow associated with 
Hq~ and Ho and the density held of figure 2: 

• the maximum of Ho takes the value of ~ 73km/(s Mpc), roughly in the same location 
as the density minimum in the void centre, denoting the fastest expansion rate in the 
central void. The maxima of Hq l * (not displayed) occur also in the central void and 
have similar magnitudes. 

• the minima of Ho in the left panel of figure 3 roughly coincide with the density maxima 
in the right panel of figure 2, in agreement with the expected slower expansion rate of ex¬ 
panding over-densities. Regarding the expansion eigenvalues at the over-density (right 
panel), the minima of Hq ' = Hq y remain positive with values close to the background 
Ho ~ 60km/(s Mpc), but the minimum of Hq 1 * becomes negative (~ —40km/(s Mpc)). 
In fact, Hq 1 * < 0 holds for all the over-densities, thus indicating unequivocally that 
these structures have started undergoing a pancake collapse at t = to- Further, as 
shown in figure 5, these structures end up collapsing into a pancake shaped shell cross¬ 
ing singularity at times much later than to- However, the virialisation process occurs 
before these singularities are approached, indicating how the description of structure 
formation by means of Szekeres models breaks down (as with the spherical collapse 
model). 

The inhomogeneous Hubble scalar in the left panel of figure 3 reveals an average deviation 
of about ~ ±10 — 20% from the background CMB based value Ho ~ 68km/(sMpc), which 

5 The eigenvalues of and h% were computed for the metric (2.1)—(2.3) . Since they are coordinate 

independent invariant quantities, they are independent of the choice of metric components. 
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0.0 0.5 1.0 1.5 2.0 

t/to 


Figure 5. Scale factors and pancake collapse. The figure displays the scale factors given by 
(6.3), plotted as as functions of t/to for the spatial coordinates of the over-density of the right panel 
of figure 3. Notice the pancake collapse of the structure in which the scale factor (associated with 
H*- 1 ') is already decreasing at t/to = 1 > while £^ = £( 3 ) (associated with iT 2 '* = H* 3 ' 1 ) are increasing 
for all t. Eventually a shell crossing singularity occurs as —> 0 at some t = £ sx ~ 3to, when the 

Szekeres description of structure formation breaks down. 


is compatible with the results of [30] in which inhomogeneities were modelled by Newtonian 
numerical simulations. 


5.3 Radial peculiar velocities. 

The radial peculiar velocities of the structures relative to the background Hubble flow (iden¬ 
tified with the CMB frame) can be computed from v p™ b = [Ho ao — %o“o](x ~ Xb), where 
Xb X* i s a sufficiently large value of the normalised comoving radius x that can be iden¬ 
tified with the asymptotic ACDM background, so that Up™ b (x;,) ~ 0. These velocities are 
depicted in the left panel of figure 4, showing low density regions expanding away from the 
background frame at v ~ — 1200 km/s, while the over-densities fall into this frame at 
~1000-1200 km/s. The latter are not comparable to our CMB dipole velocity ~ 370 km/s, 
as they are infall velocities of 10 — 20 Mpc structures into the ACDM background and thus do 
not take into account infall velocities of observers inside these structures with respect to their 
centres of mass. Instead, the peculiar velocities of the over-densities in figure 4 should be 
compared with the estimated infall velocity Wp™ b ~600 km/s of our local group with respect 
to the ACDM background. In fact, these velocities are roughly compatible with similar veloc¬ 
ities reported for the range of length scales under consideration: from data and observations 
[31, 32] and from numerical simulations [33-35]. 

Radial peculiar velocities with respect to an observer at the void centre comoving with 
the origin (depicted by the right panel) are computed from u™/ 1 = [H — 77| r= o) clqt. These 
velocities exhibit an expansion away that is roughly linearly proportional to the radial area 
distance to the void centre, reaching ~ 2000 km/s for structures located ~ 30 Mpc away. 
As shown in the right panel of figure 4, these velocities closely match the peculiar velocities 
observed [36] for galaxies in the Virgo supercluster with respect to an observer in the centre 
of the local void. 
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6 Theoretical issues 


6.1 An exact relativistic analogue of the Zeldovich approximation. 

The pancake collapse of the Szekeres structures suggests a non-trivial connection to the 
Newtonian Zeldovich Approximation (ZA) [37, 38]. Since the Szekeres models are an exact 
solution of Einstein’s equations, this issue has been examined in various attempts to obtain 
relativistic generalisation of the ZA [39-44], We provide here a brief qualitative discussion 
that reveals how the Szekeres pancake collapse yields an exact relativistic analogue of the ZA 
(for comprehensive treatment see [39-44]). 

In Newtonian gravity the “displacement” relating Eulerian y n , n = 1,2,3 and Lagrangian 
coordinates x n for a homogeneously expanding medium is y n = a(t)x n 6 . The ZA considers 
for dust sources a first order correction through a displacement y n = a(t)[x n + T n (t, x m )\ 
(with (t in i, x m ) = 0). Assuming ty n m to be symmetric we obtain the following form of the 
density [37, 38] 

_ Pini _ _ Pini _ , 

p -detyn m - { J 

where —^ A \t,x m ) (A = 1,2,3) are the eigenvalues of the “deformation” tensor = T™ 
whose components take a simple form in the coordinates x A that diagonalise ^ n m - Since for 
arbitrary displacements the three eigenvalues are different (say 0 < < ^W), 

eventually as they grow in time we have —> 1 while < 1, producing a “wall” or 

“pancake” shape deformation in which distances contract in the direction of the eigenvector 
of and expand in the directions of the other two eigenvectors, leading in the end to a 
density caustic or singularity (/? — > oo as det y n rn —» 0). 

A direct qualitative comparison of the above mentioned process can be established with 
the evolution of Szekeres models through the exact density form 


_ Pi ni _ Jini ^ _ /XT- _ a 3 (T - W) r sind 

P ~ £(1)*(2) £ (3) “ Pi J ’ J ~ Vd6t 9mn ^ _ £^2 ’ 


( 6 . 2 ) 


where g mn is the spatial metric in (2.1)-(2.2) and the normalised scale factors follow 
from the eigenvalues H ( ‘ 4 ^ = /£^ of the expansion tensor T-L^ in (5.1) with the condition 

= 1 

*Tm 

/ 1 )= a g= a j r _~^ | ^)=l^)= a . (6.3) 

Comparison of (6.1) and (6.2) yields immediately 


= 1 


£W _ a T - W 
IT ~ ”al-W’ 


^( 2 ) = ^( 3 ) = 1 


£W 

a 



(6.4) 


so that the 3-dimensional deformation matrix takes the from [S A 5^ + W]. 

where x A are the coordinates that diagonalise the spatial metric (see [17] and Appendix A 
of [23]). The Szekeres configurations we have studied provide an exact relativistic analogue 
of the ZA, as the pancake collapse of the structures occurs in an analogue manner as de¬ 
scribed before: at t = fi n i we have zero deformation = 0, but as shown in figure 5, 
we have along the over-densities 1^ > £^ as these scale factors grow for t > t\„i. so that 

°These coordinates bear no relation with the [x, y] coordinates used in figures 2, 3 and 4. 


- 12 











$A) —> 0 and p — y oo (shell crossing) may occur as T — W —> 0 at some t sx > to > ii n ; 
with Isx = ^sx' = a{t sx, r) > 0. Hence, the over-densities we have studied exhibit the same 
type of pancake deformation and collapse as Newtonian structures studied by means of the 
ZA, leading also to a final singularity (shell crossing) as —> 1 occurs as t —>• t sx , while 

0 < £sx = dx < 1 holds. 


6.2 Connection with linear perturbations. 

When considering the early evolution of inhomogeneities, the deviation from an FLRW back¬ 
ground is small and we may present the inhomogeneities as perturbations of the otherwise 
ACDM homogeneous universe, itself described by the background quantities A(t) = {p, 77, 1C}. 
The perturbative description is valid for a regime where, given a small positive parameter 
e <C 1, the following relations hold at a given initial time time fi n ; = t LS : 

\A q (r,t hS ) - A(t LS )| « 0(e), \rA' q (r, t LS )| « 0(e). (6.5) 

Since at the initial time we assume a LS = T LS = 1, the above conditions imposed on Eqs. (2.14) 
and (2.15) imply, 

~ e, |a[s ) (a 6 i ,^)| « e \A hS (r,0,(j>) - A hS \ & e, (6.6) 


Further, it can be shown that for cosmic time intervals sufficiently close to t LS the metric 
variables a and T will satisfy for all r that (see proof in [27]) 

a — a(t ) ~ 0(e), F — 1 ~ 0(e). (6-7) 

As a consequence of (6.5)— (6. 7) , up to first order in e the evolution equations (2.5) and 
(2.6) are identical to the energy conservation and Raychaudhuri equations of the ACDM 
background, while the constraint (2.11) is the background Friedman equation and (2.9)- 
(2.10) are the definitions of 77 and D^. The remaining evolution equations (2.7) and (2.8) 
and the constraint (2.12) can then be linearised up to first order in e, leading to: 


A (p) = -3D (w) + 0(e 2 ) 

(6.8) 

D (n) = -277 D (w) - —pA (p) + 0(e 2 ), 

3 

(6.9) 

D (k) = —pA (p) - 277D (h) + 0(e 2 ), 

3 

(6.10) 

which combine into the second order equation: 

AW + 277 AW - 4vrpAW + 0(e 2 ) = 0. 

(6.11) 


Evidently, (6.8)-(6.11) are mathematically equivalent to the linear dynamical equations of the 
Cosmological Perturbation Theory (CPT) formalism in the isochronous gauge (cf. [45, 46]). 
However, the exact fluctuations A ^ and relate to the gradients of A q via (2.14)—(2.15), 
and thus are analogous but not strictly equivalent to CPT perturbations. The rigorous 
equivalence between Szekeres scalars and CPT variables follows from extending to Szekeres 
the results obtained in [27] for LTB models. Instead of the q-scalars A q in (2.13), we consider 
the functional averages (A) q [?’;,] in bounded comoving domains V[rb\ along the time slices, 
which define the non-local exact fluctuations 


A 


(p) 

NL 


M 


p{t,x J ) - ( p) q [rb](t ) 

(p)g[n] 


D $ =A(t,xi)-(A)[r b ](t), A = p,H,JC , (6.12) 
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For models admitting an asymptotic ACDM background, the scalars of the CPT background 
become rigorously defined by the functional averages (A) q [r^\ evaluated in asymptotic domains 
extending to the complete time slices through the asymptotic limit 

lim {{p)q[r b \, (%)[r 6 ], {1C) in]} = {Pas,^as,^as} = { p, U , ic } . (6.13) 


while the non-local exact fluctuations in (6.12) lead in this limit to the asymptotic exact 
fluctuations: 


*(p) = ~ Pas(t) 

p as (f) 


lim 

r^-too 


M, 


= A(t, x 3 ) — A aa (t) 


lim 

rt,—>oo 


In] 

(6.14) 


so that = 5 = p/p — 1 (but ^ A^), where 5 is the density contrast plotted in figure 
2 and obtained from the exact Szekeres density 1 + 6 = [pq/p^] [1 + A^j. Actually, the 
conditions in eqs. (6.6) and (6.7) guarantee that the difference is of order e and consequently, 
their amplitude too: 


Ah’) « A (p) , D (w) « D (w) , D (k:) « D (/c) . (6.15) 

The evolution equations for these variables are mathematically identical to (6.8)-(6.11) up 
to order e (see [27]). As a consequence, the evolution of the linear order quantities is exactly 
the well known evolution of the dust perturbations in the synchronous and comoving gauge 
(see e.g. [47]). 

The equivalences presented in this subsection show unequivocally that the exact inho¬ 
mogeneities of the Szekeres models can be directly connected with the perturbative approach 
to the study of large scale structure formation. 

7 Conclusions. 

We have shown how the dynamical freedom of Szekeres models makes it possible to obtain 
a fully relativistic non-perturbative description of non-trivial networks of cold dark matter 
structures (over-densities and density voids) evolving from the last scattering surface to the 
present. In particular, we provided a numerical example of the evolution of a ~ 80Mpc sized 
region immersed in a ACDM background, consisting of a spheroidal density void surrounded 
by ten pancake shaped density maxima, placed at given radial and angular comoving locations, 
specified by the initial conditions displayed in Table 1. This configuration (whose density 
contrast is depicted in figure 2) represents a huge improvement over previous attempts to 
model cosmic structure with LTB models [4, 5, 48] or Szkeres models of class I [6, 7, 16] and 
class II [49], as they furnish a significantly better (though still coarse grained) description 
of cosmic structures observed, or inferred, at a ~100 Mpc scale today (as for example in 
[24, 25]). 

By looking at the Hubble scalar H and the eigenvalues of the expansion tensor H^ 1 ', H (2 ^ = 
H®, we have shown through a numerical example how an ~ 80 Mpc region that expands on 
average (7~L > 0 for all t) contains various ~ 10-20 Mpc over-densities undergoing a local “pan¬ 
cake” collapse at present cosmic time (H, 1 , 1 '* < 0, h< 2) = h< 3 > > 0). We have also examined 
for these structures the radial peculiar velocities with respect to the CMB background frame 
and with respect to an observer in the void centre. These velocities fit very well observed 
velocities reported in the literature for same scale structures. 
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We have also shown how the resulting Szekeres models relate to other standard theo¬ 
retical frameworks considered in current cosmological research. Specifically, we show (i) how 
the pancake collapse of the Szekeres over-densities provides an exact relativistic analogue to 
the pancake collapse in the Newtonian Zeldovich approximation, and (ii) how their evolution 
in the linear regime relates to cosmological perturbations of dust sources in the synchronous 
gauge. 

As mentioned before (last paragraph of section 3), our models are not related to Swiss 
Cheese models. Rather, they generalise the “walls and voids” models examined in [50] and the 
“onion” models of [51], all based on LTB solutions (though the onion model has been applied to 
simple dipole Szekeres models in [21, 22]). Since the scalar averaging of non-spherical Szekeres 
scalars is spherically symmetric [17] (in Buchert’s formalism and in quasi-local averaging), 
then the scalar averaging of the structures we have studied should yield similar results as 
the averaged LTB walls and voids models of [50]. However, this issue needs to be carefully 
verified and thus will be examined in a separate article. 

While the possibility of describing the relativistic and non-perturbative evolution of 
elaborated networks of non-spherical over-densities and voids is very appealing, the Szekeres 
models we have studied exhibit the expected limitations characteristic of all analytic or semi- 
analytic structure formation models. Evidently, it is wholly unreasonable to expect these 
models to describe the virialisation of cosmic structures or complex dynamical interactions 
like the “bullet cluster” or mergers of structures. As the over-densities undergo a pancake 
collapse the models break down when the shell crossing singularity is approached, and thus 
the description of the virialisation process must be introduced “by hand” as is done with 
the spherical collapse model. Still, notwithstanding the above mentioned limitations, these 
models have an enormous potential for application in open problems of current cosmological 
research, such as: 

1. Exploring the effects of relativistic corrections in cosmographic studies [24, 25] and 
Newtonian simulations (e.g. [26]), as well as the effects of non-linearity in perturbative 
relativistic treatments [1-3]. 

2. Addressing the apparent tension between estimations of the Hubble constant from the 
CMB and from supernovae surveys [52], and the associated problem of the differential 
expansion produced by nearby nonlinear structures [53-55]. 

3. Verifying if Szekeres models allow for the description of collapse regimes besides pan¬ 
cake collapse (spherical or filamentary). Explore in more detail the connection between 
evolving Szekeres structures and fully relativistic generalisations of the Zeldovich ap¬ 
proximation [39-44, 57, 58] and implementing them [1-3] in specific structure formation 
scenarios. 

4. A nonlinear relativistic treatment of peculiar velocities [56-58] and interpretation of 
the observed redshift space distortions [59, 60]. This is absent in the literature and the 
Szekeres models can help to fill this gap too. 

5. Exploring the relativistic corrections in structure formation scenarios examined by 
means of Newtonian simulations [61-63] and those attributed to modified gravity [64, 
65] , as well as the correspondence and equivalence of exact solutions vs linear and non¬ 
linear perturbative approaches [66, 67] or an extension of previous work on these issues 
in LTB models [27]. 
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6. The Szekeres models we have examined can serve as non-trivial exact “test models” 
to probe important issues theoretical, such as back-reaction, averaging and the “fitting 
problem” [68-71]. Besides looking at a more realistic non-spherical extension of the 
results of [50], we can use the Szekeres configurations as “test models” to examine 
controversial theoretical issues on back-reaction [72-75]. 

These possible applications will be pursued in separate articles currently under elaboration 
for future submission. 
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